Hidden Markov Models: Basics for animal movement data
Marie Auger-Méthé & Ron Togunov
1 Tutorial goals and set up
The goal of this tutorial is to explore how to fit hidden Markov models
(HMMs) to movement data. To do so, we will investigate a new R package,
momentuHMM. This package builds on a slightly older package,
moveHMM, that was developed by Théo Michelot , Roland Langrock, and
Toby Patterson, see associated paper:
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12578.
momentuHMM, was developed by Brett McClintock and Théo Michelot.
momentuHMM has new features such as allowing for more data streams,
inclusion of covariates as raster layers, and much more, see associated
paper:
https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.12995.
The primary learning objectives in this first tutorial are to:
- Fit simple HMMs using
momentuHMM - Checking model fit
- Determining number of states to use
- Incorporating and interpreting covariates on behaviour transition probabilities
- Exploring effect of using different temporal resolutions
1.1 Setup and data preparation
First, let’s load the packages that we will need to complete the analyses. Off course you need to have them installed first.
library(momentuHMM) # Package for fitting HMMs, builds on moveHMM
library(raster) # For importing and extracting raster spatial covariatesWhile the raster package is being phased out in favour of the new terra package, momentuHMM still relies on the raster package to extract covariates. Therefore, for this tutorial we will still use raster, but do start working with terra in your work where possible for future compatability.
Make sure to set working directory to “Day1” of the HMM workshop folder:
One of the main features of the data used with HMMs is that locations are taken at regular time steps and that there are negligible position error. For example, HMMs are adequate to use on GPS tags that take locations at a set temporal frequency (e.g., every 2 hrs). Without reprocessing, HMMs are not particularly good for irregular times series of locations or locations with large measurement error (e.g., you would need to preprocess Argos data before applying HMMs).
Unfortunately, movement of aquatic species is rarely taken at regular time intervals and without measurement error. For example, the data set we will work on, is the movement track of a grey seal tagged with a Fastlock GPS tag. This is a subset of the data set used in the paper from Whoriskey et al. 2017 (https://onlinelibrary.wiley.com/doi/abs/10.1002/ece3.2795), which applies a set of different HMMs to the movement data of various aquatic species. The original dataset is available in the online supplementary information.
Let’s read the data file and peak at it.
# Make sure you are in the good directory!
seal <- read.csv("data/seals.csv", stringsAsFactors = FALSE)
# Look at the data
head(seal)## date lon lat
## 1 2014-02-18 03:32:08 -60.16754 43.97577
## 2 2014-02-18 15:32:08 -60.16785 43.97579
## 3 2014-02-19 03:32:08 -60.16801 43.97568
## 4 2014-02-19 11:00:07 -60.17035 43.97532
## 5 2014-02-19 12:02:40 -60.20350 43.94928
## 6 2014-02-19 13:04:07 -60.22503 43.91563The files has three columns: date with the date and time, lon with the
longitude, and lat with latitude. It’s clear from the date column that
the locations are taken at irregular time intervals. So just like in
Whoriskey et al. 2017, we will regularize it. Note, that momentuHMM
has functions to preprocess irregular movement time series and/or time
series with large measurement errors (e.g., crawlWrap).
However, these
functions are much more complex.
We will cover some of these functions in tomorrow’s tutorial.
For now, as in
Whoriskey et al. 2017, we will assume that the animal moves in a
straight line and at constant speed between two locations and
interpolate the location at regular time intervals based this
assumption. Note that this kind of interpolation may not be adequate in
many circumstances and see activities to investigate the effects of the
choice of the time interval on the analysis.
# First, let's transform the date/time info into a
# proper time format
seal$date <- as.POSIXct(seal$date, format = "%Y-%m-%d %H:%M:%S",
tz = "GMT")
# Let's find the most appropriate time interval
# Calculate the time intervals
tint <- as.numeric(diff(seal$date), units = "hours")
# Let's look at the time intervals
summary(tint)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.9997 1.0131 1.0579 1.6664 2.0333 25.0017# Min. 1st Qu. Median Mean 3rd Qu. Max.
plot(tint, ylab = "Time interval (hr)", pch = 19, cex = 0.5,
las = 1)
# Let's go for 2 hrs, look at the activities to see
# the effect of this choice
# Regularising Create a 2hr sequence (2*60*60) from
# the first time to the last time
ti <- seq(seal$date[1], seal$date[length(seal$date)],
by = 2 * 60 * 60)
head(ti)## [1] "2014-02-18 03:32:08 GMT" "2014-02-18 05:32:08 GMT"
## [3] "2014-02-18 07:32:08 GMT" "2014-02-18 09:32:08 GMT"
## [5] "2014-02-18 11:32:08 GMT" "2014-02-18 13:32:08 GMT"# Interpolate the location at the times from the sequence
iLoc <- as.data.frame(cbind(lon = approx(seal$date, seal$lon,
xout = ti)$y, lat = approx(seal$date, seal$lat, xout = ti)$y))
# Create a new object that has the regular time and
# interpolated locations
sealreg <- cbind(date = ti, iLoc)
head(sealreg)## date lon lat
## 1 2014-02-18 03:32:08 -60.16754 43.97577
## 2 2014-02-18 05:32:08 -60.16759 43.97577
## 3 2014-02-18 07:32:08 -60.16764 43.97578
## 4 2014-02-18 09:32:08 -60.16770 43.97578
## 5 2014-02-18 11:32:08 -60.16775 43.97578
## 6 2014-02-18 13:32:08 -60.16780 43.97578# Quickly plot the regularized locations
plot(sealreg$lon, sealreg$lat, pch = 19, cex = 0.5, xlab = "Lon",
ylab = "Lat", las = 1)
# Add the original data
points(seal$lon, seal$lat, pch = 19, cex = 0.5, col = rgb(1,
0, 0, 0.5))
# add legend
legend(x = -62.4, y = 42.1, pch = c(19, 19),
legend = c("Original data", "Regularised locations"),
col = c("red", "black"))
Here, we interpolated the location data to 2 hours. The resolution we choose can have a significant effect on the results, so the this decision must be made wisely. At the end of this tutorial, in section Effects of interpolation. In tomorrow’s tutorial, we will review a more comprehensive strategy to select an appropriate temporal resolution when we have irregular data.
Ok, we have now a regular time series of locations. For the most basic
movement HMMs, you need to transform the time series of locations into
two time-series: one which represents the step lengths (distance between
two locations) and one that represents the turning angle (angle between
two steps). The package momentuHMM has a function that does that:
prepData. The data is in latitude and longitude and I’m assuming using
WGS84. You can use prepData on latitude and longitude data or on
projected data (e.g., UTM). prepData can do much more, and we will see
some other features later, but for now the arguments that we need to
think about are:
type: whether it’s easting/northing coordinate system (
type="UTM") or whether it’s longitude/latitude (type = "LL"). We are using"LL"."coordNames": the names of the columns with the coordinates. So in our caselonandlat.
# Preparing the data: here mainly calculating step
# lengths and turning angles
sealPrep <- prepData(sealreg, type = "LL", coordNames = c("lon",
"lat"))
# Let's peak at the data
head(sealPrep)## ID step angle date x y
## 1 Animal1 0.004147236 NA 2014-02-18 03:32:08 -60.16754 43.97577
## 2 Animal1 0.004147236 -4.615664e-09 2014-02-18 05:32:08 -60.16759 43.97577
## 3 Animal1 0.004147236 -4.567300e-09 2014-02-18 07:32:08 -60.16764 43.97578
## 4 Animal1 0.004147235 -4.124921e-09 2014-02-18 09:32:08 -60.16770 43.97578
## 5 Animal1 0.004147235 -5.074839e-09 2014-02-18 11:32:08 -60.16775 43.97578
## 6 Animal1 0.004147235 -4.346000e-09 2014-02-18 13:32:08 -60.16780 43.97578## [1] "momentuHMMData" "data.frame"If we take a quick look at the data, we can see that the step lengths
(column step) and angles (column angle) have been calculated and that
the names of the columns with the coordinates are now x and y, rather
than lon and lat. We also have a new column: ID. In cases where
you have multiple individuals, you would want your original data to have
an ID column before being imputed in the prepData function.
Note that the object it returns is from a specific S3 class defined in
the momentuHMM package. This means that we can apply generic functions
like plot to it and it will return specified outputs.
We can see the track of the animal and the times series of the step lengths and turning angles, as well as the histograms of the step lengths and turning angles. Can you already identify patterns?
1.2 Fitting the HMM
Now that our data is ready to use, we can fit the HMM. For this we use
the function fitHMM. This is where we need to make many of our
modelling decisions and most of these decisions will be associated with
one of the arguments of the function. The minimum arguments fitHMM
requires are: fitHMM(data, nbState, dist, Par0).
First, when we fit a HMM, we need to decide the number of behavioural
states we are going to model. To start simple, we will only use two
behavioural states. These could be, for example, representing one
behaviour with fast and directed movement (e.g., travelling) and one
behaviour with a more slow and tortuous movement (e.g., foraging). This
mean that the argument nbState will be set to 2. Second, we need to
decide whether we make the transition probabilities between the
behavioural states dependent on covariates. Here again, we will start
simple and we will not use covariates. That means that our formula
argument will be ~1. Third, we need to select the distribution we will
use to model the step lengths and the one we will use to model the
turning angles. For now, we will use the gamma distribution for the step
length and the von Mises for the turning angles. This means that the
argument dist will be set to: list(step="gamma", angle="vm"). Note
that dist should be a list with an item for each data stream columns
in the data that we are modelling (so here the column step and
angle). The gamma distribution is strictly positive (i.e., it does not
allow for 0s). If you have step lengths that are exactly zero in your
data set, you need to use zero-inflated distributions. But in this case,
we have no zeros. Let’s check
## [1] 0Since the sum is 0, we don’t have any values of zero and don’t need to worry about modelling zero-inflation.
By default, fitHMM will set the argument estAngleMean to NULL,
which means that we assume that the mean angle is 0 for both behaviours
(i.e., the animal has a tendency to continue in the same direction) and
that only the angle concentration parameters differ. The concentration
parameters control the extent to which the animal continues forward
versus turn. Doing so reduces the number of parameters to estimate.
These are all very important decisions that you must make when you
construct your model. In addition, we need to specify initial values for
the parameters to estimate. The HMMs are fit using maximum likelihood
estimate (MLE) with a numerical optimizer. An unfortunate aspect of
fitting models using numerical optimizers, is that, to be able to
explore the parameter space and find the best parameter values for the
model, the optimizer needs to start somewhere. You need to decide where
it starts. Unfortunately, choosing bad starting values can result in
parameter estimates that are not the MLE, but just a local maximum. To
choose your initial parameter you can take a quick peak at the data
(e.g., using the plots below) and use some general biological
information. For example, it’s common for animals to have one behaviour
with long step lengths and small turning angles and one behaviour with
short step lengths and larger turning angles. From the plots above, it
looks like the animal has step lengths that are close to 0, 3, 7. As I
see it, there are probably three behaviours with different mean step
lengths around these values.
plot(sealPrep$step ~
sealPrep$date, ty = "l", ylab = "Step length",
xlab = "Date", las = 1)
abline(h = 0, col = rgb(1, 0, 0, 0.5))
abline(h = 3, col = rgb(1, 0, 0, 0.5))
abline(h = 7, col = rgb(1, 0, 0, 0.5))
So let’s choose 3 and 7 for the means step lengths and use the same values for the standard deviations. The turning angles are either very close to 0 or pretty spread from \(-\pi\) to \(\pi\). High concentration parameter values (\(\kappa\), said kappa) mean that the animal has a tendency to move in the same direction. Values close to 0 mean that the turning angle distribution is almost uniform (the animal turns in all directions). Note that \(\kappa\) cannot be exactly 0, so let’s choose 0.1 and 1.
# Setting up the starting values
mu0 <- c(3, 7) # Mean step length
sigma0 <- c(3, 7) # Sd of the step length
kappa0 <- c(0.1, 1) # Turning angle concentration parameter (kappa > 0)Ok, were are ready. Let’s fit the HMM
# Fit a 2 state HMM
sealHMM2s <- fitHMM(sealPrep, nbState = 2, dist = list(step = "gamma",
angle = "vm"), Par0 = list(step = c(mu0, sigma0), angle = kappa0), formula = ~ 1)Let’s explore the results.
## Value of the maximum log-likelihood: -1810.643
##
##
## step parameters:
## ----------------
## state 1 state 2
## mean 0.05014675 4.902540
## sd 0.06428477 2.689235
##
## angle parameters:
## -----------------
## state 1 state 2
## mean 0.000000 0.000000
## concentration 2.789104 1.771595
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.346325 -4.840844
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2
## state 1 0.91264169 0.08735831
## state 2 0.00783846 0.99216154
##
## Initial distribution:
## ---------------------
## state 1 state 2
## 9.999995e-01 5.159521e-07## Decoding state sequence... DONE
Based on the mean step parameters, it looks like the first behavioural state (state 1) has really small step lengths compare to state 2. This is particularly easy to see in the step histogram (first figure), where the estimated distribution for each state is overlaid on top of the observed step length frequencies. Surprisingly, the turning angle distributions of the two states both indicate directed movement (second figure). In fact, the concentration parameter of state 1 is bigger than that of state 2. In the last figure with the track, it looks like if we only have state 2. This is strange, since the transition probabilities estimated indicate that the animal remain in each behaviour for multiple steps (diagonal values of the transition probability matrix are > 0.5).
Part of what’s happening may be that the locations of state 1 are hidden under state 2. We can see if this is the case by plotting state 2 with a higher transparency.
# plot it with transparency
plot(sealHMM2s, col = c(rgb(230, 160, 20, max = 255),
rgb(80, 180, 230, alpha = 80, max = 255)))## Decoding state sequence... DONE
Indeed, we are getting two states, but there is likely more going on.
2 Identifying behavioural states
Maybe looking at the state sequence in more detail will help us
understand. Actually, we are interested in identifying when an animal is
in each of the behavioural states (here when the animal is in state 1 vs
state 2), something we call state decoding. To do so, we can use the
function viterbi, which uses the Viterbi algorithm to produce the most
likely sequence of states according to your fit model and data.
# Apply the Viterbi algorithm using your fited model object
sealStates <- viterbi(sealHMM2s)
# Let's look at predicted states of the first 20 time steps
head(sealStates, 20)## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2## sealStates
## 1 2
## 55 504So they are some state 1, we are just not seeing them!
In many cases, it is more interesting to get the probability of being in
each state rather than the most likely state. To do so, you can use the
function stateProbs, which returns a matrix of the probability of
being in each state for each time step.
# Calculate the probability of being in each state
sealProbs <- stateProbs(sealHMM2s)
# Let's look at the state probability matrix
head(sealProbs)## state 1 state 2
## [1,] 1 9.447481e-12
## [2,] 1 7.208529e-12
## [3,] 1 7.208528e-12
## [4,] 1 7.208527e-12
## [5,] 1 7.208526e-12
## [6,] 1 7.208525e-12We can see here the probability of being in both states for the first 6 time steps. Here, the probability of being in state 1 is really high for these steps, but that may not always be the case. Sometimes you might have values closer to 0.5, which would indicate that for that time step, it’s hard to identify which state you are in (i.e., which step length and turning angle distributions fit best).
You can plot the results from both of these functions using the function
plotState
## Decoding state sequence... DONE
## Computing state probabilities... DONE
We clearly have two behavioural states. The steps from state 1 are
just too small for us to see them in the track.
3 But do we have a good model?
This is all great, but is this a good model for our data?
3.1 Starting values: should we care?
The first thing we need to look into is whether the parameter estimates are reliable. As I mentioned above, the initial parameter values can affect the estimating procedures. So it’s always good to check if you get similar results with different starting values. Here, we will make the starting values of the two behavioural states closer to one another.
# Setting up the starting values (more similar values)
mu02 <- c(5, 7) # Mean step length
sigma02 <- c(5, 7) # Sd of the step length
kappa02 <- c(1, 1) # Turning angle concentration parameter (kappa > 0)
# Fit the same 2 state HMM
sealHMM2s2 <- fitHMM(sealPrep,
nbState = 2,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu02, sigma02), angle=kappa02))Let’s compare the two results. First let’s look at which of the two has the lowest negative log likelihood (equivalent of highest log likelihood, so closer to the real MLE). Let’s look also at the parameter estimates they each returned.
## original new
## 1810.643 1810.643## state 1 state 2 state 1 state 2
## mean 0.05014675 4.902540 0.05014668 4.902536
## sd 0.06428477 2.689235 0.06428473 2.689231## state 1 state 2 state 1 state 2
## mean 0.000000 0.000000 0.000000 0.000000
## concentration 2.789104 1.771595 2.789099 1.771601## state 1 state 2 state 1 state 2
## state 1 0.91264169 0.08735831 0.912643590 0.08735641
## state 2 0.00783846 0.99216154 0.007838454 0.99216155Looks like they both returned very close values for everything. So
that’s good! The function fitHMM also has the argument retryFits
which perturbs the parameter estimates and retry fitting the model. The
argument is used to indicate the number of times you want perturb the
parameters and retry fitting the model (you can choose the size of the
perturbation by setting the argument retrySD).
Let’s try (this will take a few minutes).
# Fit the same 2-state HMM with retryFits
# This is a random pertubation, so setting the seed to get the same results
set.seed(1234)
sealHMM2sRF <- fitHMM(sealPrep,
nbState = 2,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu02, sigma02), angle=kappa02),
retryFits=10)## Attempting to improve fit using random perturbation. Press 'esc' to force exit from 'fitHMM'
##
Attempt 1 of 10 -- current log-likelihood value: -1810.643
Attempt 2 of 10 -- current log-likelihood value: -1810.643
Attempt 3 of 10 -- current log-likelihood value: -1810.643
Attempt 4 of 10 -- current log-likelihood value: -1810.643
Attempt 5 of 10 -- current log-likelihood value: -1810.643
Attempt 6 of 10 -- current log-likelihood value: -1810.643
Attempt 7 of 10 -- current log-likelihood value: -1810.643
Attempt 8 of 10 -- current log-likelihood value: -1810.643
Attempt 9 of 10 -- current log-likelihood value: -1810.643
Attempt 10 of 10 -- current log-likelihood value: -1810.643Let’s compare the results again.
# Negative log likelihood
c(original = sealHMM2s$mod$minimum, new = sealHMM2s2$mod$minimum,
retryFits = sealHMM2sRF$mod$minimum)## original new retryFits
## 1810.643 1810.643 1810.643## state 1 state 2 state 1 state 2 state 1 state 2
## mean 0.05014675 4.902540 0.05014668 4.902536 0.05014666 4.902544
## sd 0.06428477 2.689235 0.06428473 2.689231 0.06428477 2.689237## state 1 state 2 state 1 state 2 state 1 state 2
## mean 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## concentration 2.789104 1.771595 2.789099 1.771601 2.789111 1.771598## state 1 state 2 state 1 state 2 state 1 state 2
## state 1 0.91264169 0.08735831 0.912643590 0.08735641 0.912642638 0.08735736
## state 2 0.00783846 0.99216154 0.007838454 0.99216155 0.007838521 0.99216148Still looking very similar!
But let’s choose some wild starting values to see the type of problems we might run into.
# Setting wild starting values
mu0W <- c(5, 100) # Mean step length
sigma0W <- c(5, 0.1) # Sd of the step length
kappa0W <- c(10, 0.01) # Turning angle concentration parameter (kappa > 0)
sealHMM2sW <- fitHMM(sealPrep, nbState = 2,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu0W,sigma0W), angle=kappa0W))Let’s compare the negative log likelihood (note that we want the minimum value here to get the MLE) and parameter estimates.
## original wild
## 1810.643 2122.390## state 1 state 2 state 1 state 2
## mean 0.05014675 4.902540 4.420324 100.0
## sd 0.06428477 2.689235 4.479956 0.1## state 1 state 2 state 1 state 2
## mean 0.000000 0.000000 0.000000 0.00
## concentration 2.789104 1.771595 1.840826 0.01## state 1 state 2 state 1 state 2
## state 1 0.91264169 0.08735831 1.0000000 0.000000e+00
## state 2 0.00783846 0.99216154 0.9999985 1.485305e-06The negative log likelihood here is much larger than the negative log
likelihood of the model fit with the sensible starting values. Note
also that some of the parameter estimates are exactly the same values as
the initial values, which is always a sign of optimization problems.
These are all signs that the model fit with this new set of starting
values is not good. Sometimes, you will get error messages saying that
nlm or optim, the functions that minimize the negative log
likelihood, did not converged or that non-finite values were supplied to
them. Sometimes, you might not be able to plot the data or get the
states. These are also signs that the starting values are poor or that
the model is not good for your data.
We will see later, that the package momentuHMM has functions to help
selecting initial parameters for complex models. However, these
functions rely on the user choosing initial parameter values for a
simple model like the one we just explored.
3.2 Pseudo-residuals
Ok, we fit a model to data, looks like the parameter estimates are
reliable, and we can get state probabilities, but is this a good model
for our data? Can we use the model results? The best way to investigate
model fit is through pseudo-residuals. Pseudo-residuals are a type of
model residuals that account for the interdependence of observations.
They are calculated for each of the time series (e.g., you will have
pseudo-residuals for the step length time series and for the turning
angle time series). If the fit model is appropriate, the
pseudo-residuals produced by the functions pseudoRes should follow a
standard normal distribution. You can look at pseudo-residuals directly
via the function plotPR, which plots the pseudo-residual times-series,
the qq-plots, and the autocorrelation functions (ACF).
## Computing pseudo-residuals...
Both the qq-plot and the ACF plot indicate that the model does not fit the step length time series particularly well. The ACF indicates that there is severe autocorrelation remaining in the time series, even after accounting for the persistence in the underlying behavioural states.
This could indicate that there are more hidden behavioural states. Let’s try a 3-state HMMs.
# Setting up the starting values (we need 3 now, because 3 states)
mu03s <- c(0.1, 3, 7) # Mean step length
sigma03s <- c(0.1, 3, 7) # Sd of the step length
kappa03s <- c(0.1, 1, 1) # Turning angle concentration parameter
# Fit a 3 state HMM
sealHMM3s <- fitHMM(sealPrep,
nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu03s, sigma03s), angle=kappa03s))Let’s look at it and at the pseudo-residuals.
## Decoding state sequence... DONE
## Decoding state sequence... DONE
## Computing state probabilities... DONE
## Computing pseudo-residuals...
Ah, that’s better. Not perfect, but less unexplained autocorrelation, especially in the step lengths. Some of the unexplained autocorrelation in the turning angles is related to the method we used to get location estimates at regular time intervals, see activities to explore this a bit more in depth. If we look at the step length distribution, we can see that we have still our state with really small steps (maybe resting?), a state with steps of medium size (maybe foraging?), and a state with large steps (maybe travelling?). Looking at the track, it does look like the movement in state 2 is more tortuous and in focal areas, while the movement in state 3 is very directed and span larger areas.
3.3 Simulating data
We can also use the function simData to look at whether the movement
produced by the models is similar to the observed movement. Let’s
simulate the 3-state HMMs, and then fit it to be able to identify the
behavioural states. Note that we could just plot directly the simulated
data and colour code the state, rather than estimating the fitting the model and predicting the
states and plotting the predicted states.
However, plotting the simulated data might require a bit of fiddling (A good thing to know is
that the colour palette used by the packages is:
c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")).
One way fitting the model to the simulated data may be useful is in
identifying whether you have enough data to estimate your parameters
accurately.
# Simulate the model
sealSim <- simData(model=sealHMM3s, nbAnimals =1, states = TRUE)
# Fit the model to the simulated data
sealSimFit <- fitHMM(sealSim,
nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu03s, sigma03s), angle=kappa03s))## Decoding state sequence... DONE
## Value of the maximum log-likelihood: -3474.566
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3
## mean 0.05783530 3.772843 7.352621
## sd 0.07358705 2.036328 1.455070
##
## angle parameters:
## -----------------
## state 1 state 2 state 3
## mean 0.000000 0.000000 0.000000
## concentration 2.959238 1.285337 9.212642
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.644707 -3.716304 -4.678914 -3.844702 -44.02199 -2.69548
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2 state 3
## state 1 9.129503e-01 0.06484338 0.02220633
## state 2 9.012574e-03 0.97023148 0.02075595
## state 3 7.130482e-20 0.06324057 0.93675943
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3
## 9.999993e-01 3.973350e-07 3.144418e-07## Value of the maximum log-likelihood: -1605.31
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3
## mean 0.05215803 3.809156 7.357247
## sd 0.06769963 2.018804 1.520866
##
## angle parameters:
## -----------------
## state 1 state 2 state 3
## mean 0.000000 0.000000 0.000000
## concentration 2.793543 1.244964 9.814968
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.634604 -3.758438 -4.740937 -3.638665 -4.98631 -2.777884
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2 state 3
## state 1 0.913185696 0.06551868 0.02129563
## state 2 0.008435080 0.96616689 0.02539803
## state 3 0.006389916 0.05815706 0.93545303
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3
## 9.999990e-01 4.975252e-07 4.573028e-07Looks similar, except maybe that the real seal goes back and forth to some of the same locations, something not captured in the simulation. The parameter estimates are close to those used to simulate the data, which is good.
3.4 Confidence intervals
We can look at the confidence intervals for the parameter estimates.
This does not necessarily give us an idea of whether the model is good
or not, but it helps us understand whether the behavioural states have
different movement characteristics. To compute the confidence intervals,
we can use the function CIreal.
For example, let’s look the step length mean parameters.
# Caculate the CIs
sealCI <- CIreal(sealHMM3s)
# Let's look at the step lenth CIs
rbind(lower=sealCI$step$lower["mean",],
upper = sealCI$step$upper["mean",])## state 1 state 2 state 3
## lower 0.01900231 3.572181 7.081671
## upper 0.08531376 4.046130 7.632822We can see that the 95% confidence intervals of mean step length does not overlap for the three states, indicating that, in term of step lengths, they vary quite a bit.
4 Understanding the factors that affect movement: including covariates
momentuHMM allows you to incorporate various covariates in the model
in many different ways, which allows you to understand how covariates
affect the behaviour and movement of the animal. There are so many
different options (including activity centers, spatial covariates, etc),
that it’s impossible to look at them all here. Here are a few examples.
4.1 Spatial covariates
One of the most exciting aspect of momentuHMM is that you can
incorporate spatial covariates that are kept in a raster file. Here we
will use a raster that contains the bathymetry in the area of the seal.
I created this file using the package marmap. To read the file you
need the function raster from the package with the same name.
Let’s read the file and plot it.
This represents the bathymetry in meters. Any value above 0 is on land.
Here we will use the function prepData again, but this time including
the raster as a spatial covariate. The function will extract the raster
values at the seal locations. Off course for this to work the spatial
covariate raster needs to be in the same projection as the movement
data. It’s the case here.
sealPrep <- prepData(sealreg, type="LL", coordNames = c("lon", "lat"),
spatialCovs = list(bathy=bathy))
head(sealPrep)## ID step angle date x y
## 1 Animal1 0.004147236 NA 2014-02-18 03:32:08 -60.16754 43.97577
## 2 Animal1 0.004147236 -4.615664e-09 2014-02-18 05:32:08 -60.16759 43.97577
## 3 Animal1 0.004147236 -4.567300e-09 2014-02-18 07:32:08 -60.16764 43.97578
## 4 Animal1 0.004147235 -4.124921e-09 2014-02-18 09:32:08 -60.16770 43.97578
## 5 Animal1 0.004147235 -5.074839e-09 2014-02-18 11:32:08 -60.16775 43.97578
## 6 Animal1 0.004147235 -4.346000e-09 2014-02-18 13:32:08 -60.16780 43.97578
## bathy
## 1 -18
## 2 -18
## 3 -18
## 4 -18
## 5 -18
## 6 -18Now our sealPrep object has a new column with the bathymetry value at
that location.
There are two ways in which the bathymetry can influence the movement of the animal. It can influence the behavioural state the animal is in or it can influence the movement characteristics in some of the state (e.g. change the mean step length value in one of the behavioural states).
Let’s explore the first option. So here we are assuming that the probability to switch between behavioural states is affected by the bathymetry. Because of the way the HMMs are constructed, we have to estimate one more parameter per switching probability (\(\gamma_{ij}\) said gamma i j), excluding the probability of remaining in the same state, which can be deducted based on the switching probabilities. Please see the mathy section below to understand the details. So if \(N\) represents the number of states, it means you need to estimate \(N(N-1)\) new parameters. For an HMM with 3 states, we would estimate 6 new parameters.
The main argument of the fitHMM function that we need to change are
the formula and the starting values (because we have 6 new parameters
to estimate). The default value for the formula is \(\sim 1\), which
means that the transition probability are not affected by covariates.
Here, we want to indicate that the transition probabilities are affected
by the bathymetry. We can use the function getPar0 to help us find
starting values based on our simpler 3-state HMMs. This might be a bit
of an overkill for this model, as getPar0 will set the starting values
for the covariate effects on the transition probabilities to 0. But it
can make a small difference to use good starting values for the other
parameters (might be worth looking if it makes indeed a difference).
Note that to be able to use getPar0 in this case, we have to refit the
simpler model with the new data that has the bathymetry data.
# New formula
bathyForm <- ~bathy
# Refit the simple 3-states HMM - so it's associated with the data with bathymetry
sealHMM3s <- fitHMM(sealPrep,
nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu03s, sigma03s), angle=kappa03s))
# Get new starting values based on simpler model, here the 3-states HMM
par0b <- getPar0(model=sealHMM3s, formula=bathyForm)## $Par
## $Par$step
## mean_1 mean_2 mean_3 sd_1 sd_2 sd_3
## 0.05215803 3.80915581 7.35724678 0.06769963 2.01880448 1.52086561
##
## $Par$angle
## concentration_1 concentration_2 concentration_3
## 2.793543 1.244964 9.814968
##
##
## $beta
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.634604 -3.758438 -4.740937 -3.638665 -4.98631 -2.777884
## bathy 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
##
## $delta
## [1] 9.999990e-01 4.975252e-07 4.573028e-07# Now we fit the new model with bathymetry
sealHMMb <- fitHMM(sealPrep, nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = par0b$Par, beta0=par0b$beta,
formula=bathyForm)## Value of the maximum log-likelihood: -1586.374
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3
## mean 0.05048508 3.839015 7.381269
## sd 0.06478315 2.032188 1.535275
##
## angle parameters:
## -----------------
## state 1 state 2 state 3
## mean 0.000000 0.000000 0.000000
## concentration 2.795529 1.260601 9.851826
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.20107520 -3.890246 0.5580236 -1.11891661 28.54185 -2.7256944291
## bathy 0.02380705 17.833869 0.4211101 0.02023555 28.43770 0.0009176851
##
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
## state 1 state 2 state 3
## state 1 9.970326e-01 0.002967375 0.00000000
## state 2 2.867348e-28 0.985116040 0.01488396
## state 3 0.000000e+00 0.053906987 0.94609301
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3
## 9.999978e-01 8.116212e-07 1.398467e-06The first thing to note is that the model now has a new row of
parameters in the section
Regression coeffs for the transition probabilities. This row
represents the effects of bathymetry on the transition probabilities. A
positive value means that increasing values of the covariate increase
the switching probability, while a negative value means the reverse
relationship. So for example, as bathymetry increases, the probability
of switching from state 1 to state 2 only slightly increases, but the
probability of switching from state 3 to state 1 increases much more.
Note that the relationship between these coefficients and the final
transition probabilities are dependent on each other and on the
intercept, see math section for details. Keep in mind here that
bathymetry is represented in negative values, so as bathymetry value
increases the seal is in shallower water.
The row with the intercept represents the base probabilities before we incorporate the bathymetry effects. In model without covariates, you would only get the intercept values and would use them alone to calculate the transition probabilities, see the math section below for more details.
Let’s look at it visually. First by using the generic plot function
and, then, by using the function plotStationary, which plots the
stationary state probabilities.
## Decoding state sequence... DONE
We can see that the bathymetry mainly affects the transition probability in high bathymetry values (so, close/on land). One of the most obvious patterns is that it increases the probability of switching to state 1 from both states 2 and 3 and decreases the probability of staying in state 2 and staying in state 3. Remember, state 1 is associated with the very small step lengths. This is likely seal hauling out on land or resting in shallow waters! In fact, one of the places where we have a lot of state 1 is on small islands, where I’m pretty sure they are known to haul out.
Here a quick plot to convince you.
# Create a new file, so we can see land better
water <- bathy
water[water >= 0] <- NA
# Create color ramp to show bathymetry in blue gradients
bluePal <- colorRampPalette(c(rgb(0.1,0.1,1), rgb(0.9,0.9,1)))
# plot depth with land (NA) as grey
plot(water, col = bluePal(50), colNA = "grey30")
# Let's plot only the state 1
points(sealPrep[viterbi(sealHMMb) == 1, c("x","y")],
col = rgb(230, 160, 20, 12, max = 255),
pch = 19)
# add legend
legend(x = -60.8, y = 40.3, pch=c(19, NA),
legend = c("State 1", "Land"), col = c("#E69F00", NA),
fill = c(NA, "grey20"), border = NA)
Ok, but it is a better model?
## Model AIC
## 1 sealHMMb 3218.748
## 2 sealHMM3s 3244.619Yes.
What about the pseudo-residuals?
## Computing pseudo-residuals...
Still some issues, but not too too bad.
4.2 Non spatial covariates
Could there be diurnal patterns explaining the autocorrelation in
turning angles? We are not limited to incorporate spatial covariates. We
can incorporate any covariates as long as we can relate each of the
locations to a covariate value (that’s in essence what the prepData
function does for the spatial covariates).
So here we are going to create a new column that keeps track of the hour of the day for each observation. We are going to set it in the time zone where the seals are.
Here again we can make the covariate affect the transition
probabilities. The difference is that we have to use a more complex
function that represents the cyclical nature of time of day. This is the
cosinor function. For more information on the cosinor function, see
mathy section below.
# New formula that says that the transition probabilties are affected by the time of day
diurnForm <- ~ cosinor(hour, period = 24)
# Refit model with new data
sealHMM3s <- fitHMM(sealPrep,
nbState = 3,
dist = list(step = "gamma", angle = "vm"),
Par0 = list(step = c(mu03s, sigma03s), angle = kappa03s))
# Get initial parameters
par0diurn <- getPar0(model = sealHMM3s, formula = diurnForm)
# Now we fit the diurnal model with the starting values returned by getPar0
sealHMMdiurn <- fitHMM(sealPrep,nbState = 3,
dist = list(step = "gamma", angle = "vm"),
Par0 = par0diurn$Par, beta0 = par0diurn$beta,
formula = diurnForm)## Value of the maximum log-likelihood: -1591.808
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3
## mean 0.05161135 3.859728 7.396961
## sd 0.06660846 2.043111 1.528737
##
## angle parameters:
## -----------------
## state 1 state 2 state 3
## mean 0.000000 0.000000 0.00000
## concentration 2.784346 1.269179 10.11164
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3
## (Intercept) -3.4102909 -266.1457 -4.9321125 -9.726472
## cosinorCos(hour, period = 24) 0.1381726 222.1870 -0.6235176 -5.935061
## cosinorSin(hour, period = 24) 2.1911625 152.4522 0.6596884 5.112540
## 3 -> 1 3 -> 2
## (Intercept) -252.4865 -3.0804812
## cosinorCos(hour, period = 24) 188.8899 -0.5317057
## cosinorSin(hour, period = 24) 164.8816 0.4483178
##
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
## state 1 state 2 state 3
## state 1 9.713296e-01 0.02867043 4.885871e-212
## state 2 1.306616e-02 0.96385354 2.308031e-02
## state 3 1.329916e-191 0.07286012 9.271399e-01
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3
## 1.000000e+00 3.885797e-62 9.292798e-62## Decoding state sequence... DONE
## Model AIC
## 1 sealHMMb 3218.748
## 2 sealHMMdiurn 3241.615I’m not a 100% sure how to interpret the results biologically. But again it appears to be affecting state 1, which is likely resting. According to AIC it’s less good than the bathymetry model.
To help understand the relationship further we can also use the function
plotStationary, which display the relationship between the covariate
and the stationary probability of being in each state rather than the
relationship between the covariates and the switching probability.
So it looks like there is a higher probability of being in state 2 (the foraging type of behaviour) at night and a peak in each of the other behaviours during different period of the day. Some datasets are very high resolution, and it’s important to note that daily effects can only be used when the data spans at least several days.
Covariates can affect the state-dependent movement, by which I mean the
characteristics of the movement in a given state. So for example, the
animal might make all the behaviours at night and during the day, but at
night it might do them with smaller steps. Here, we are going to look at
whether the seal is more or less directed at night. To do this, we are
going to change the argument DM. For this kind of model, it’s
particular useful to use getPar0, because it helps set starting values
for the parameters associated with the effects of covariates on the
movement parameters.
# Use the formula to describe the changes in the design matrix
diurnDM <- list(angle = list(concentration = ~ cosinor(hour, period = 24)))
# Get parameter starting values
par0diurnDM <- getPar0(model=sealHMM3s, DM=diurnDM)
# Look at starting values
par0diurnDM## $Par
## $Par$step
## mean_1 mean_2 mean_3 sd_1 sd_2 sd_3
## 0.05215803 3.80915581 7.35724678 0.06769963 2.01880448 1.52086561
##
## $Par$angle
## concentration_1:(Intercept)
## 1.0273107
## concentration_1:cosinorCos(hour, period = 24)
## 0.0000000
## concentration_1:cosinorSin(hour, period = 24)
## 0.0000000
## concentration_2:(Intercept)
## 0.2191064
## concentration_2:cosinorCos(hour, period = 24)
## 0.0000000
## concentration_2:cosinorSin(hour, period = 24)
## 0.0000000
## concentration_3:(Intercept)
## 2.2839085
## concentration_3:cosinorCos(hour, period = 24)
## 0.0000000
## concentration_3:cosinorSin(hour, period = 24)
## 0.0000000
##
##
## $beta
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.634604 -3.758438 -4.740937 -3.638665 -4.98631 -2.777884
##
## $delta
## [1] 9.999990e-01 4.975252e-07 4.573028e-07# Now we fit the new diurnal model with the starting values returned by getPar0
sealHMMdiurnDM <- fitHMM(sealPrep, nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = par0diurnDM$Par, beta0=par0diurnDM$beta,
DM=diurnDM)## Value of the maximum log-likelihood: -1594.656
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3
## mean 0.04650116 3.798213 7.362677
## sd 0.05842415 2.024625 1.514162
##
## Regression coeffs for angle parameters:
## ---------------------------------------
## concentration_1:(Intercept) concentration_1:cosinorCos(hour, period = 24)
## [1,] 1.196787 1.601231
## concentration_1:cosinorSin(hour, period = 24) concentration_2:(Intercept)
## [1,] -0.8483791 0.2166047
## concentration_2:cosinorCos(hour, period = 24)
## [1,] -0.01211124
## concentration_2:cosinorSin(hour, period = 24) concentration_3:(Intercept)
## [1,] 0.09172744 2.325207
## concentration_3:cosinorCos(hour, period = 24)
## [1,] 0.1643344
## concentration_3:cosinorSin(hour, period = 24)
## [1,] 0.2633747
##
## angle parameters (based on mean covariate values):
## --------------------------------------------------
## state 1 state 2 state 3
## concentration 0.6608244 1.258335 8.705611
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.584947 -3.907193 -4.754122 -3.582031 -4.985558 -2.740689
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2 state 3
## state 1 0.912827789 0.06882726 0.01834495
## state 2 0.008313216 0.96484561 0.02684117
## state 3 0.006380632 0.06022788 0.93339149
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3
## 9.999989e-01 5.853988e-07 5.538958e-07## Decoding state sequence... DONE
## Model AIC
## 1 sealHMMb 3218.748
## 2 sealHMMdiurnDM 3235.313
## 3 sealHMMdiurn 3241.615So here it looks like the seal makes more directed movement in state 1 at night, makes more directed movement in state 2 early in the morning, and makes more directed movement in state 3 early in the morning.
Note that we can combine various covariates. For example, here we are going to combine bathymetry and time of day.
# Get parameter starting values
par0bd <- getPar0(model=sealHMM3s, formula=bathyForm, DM=diurnDM)
# Now we fit the new model with the starting values returned by getPar0
sealHMMbd <- fitHMM(sealPrep, nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = par0bd$Par, beta0=par0bd$beta,
formula=bathyForm,
DM=diurnDM)## Value of the maximum log-likelihood: -1573.811
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3
## mean 0.04628763 3.814690 7.375478
## sd 0.05798717 2.027531 1.521130
##
## Regression coeffs for angle parameters:
## ---------------------------------------
## concentration_1:(Intercept) concentration_1:cosinorCos(hour, period = 24)
## [1,] 1.202684 1.580226
## concentration_1:cosinorSin(hour, period = 24) concentration_2:(Intercept)
## [1,] -0.8329418 0.223067
## concentration_2:cosinorCos(hour, period = 24)
## [1,] -0.009035042
## concentration_2:cosinorSin(hour, period = 24) concentration_3:(Intercept)
## [1,] 0.0889006 2.327928
## concentration_3:cosinorCos(hour, period = 24)
## [1,] 0.1680636
## concentration_3:cosinorSin(hour, period = 24)
## [1,] 0.2653602
##
## angle parameters (based on mean covariate values):
## --------------------------------------------------
## state 1 state 2 state 3
## concentration 0.6789645 1.262561 8.697037
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.47759008 -3.5164541 0.5995063 -1.68027618 20.89242 -2.699083520
## bathy 0.01086981 0.1225208 0.6682980 0.01467615 20.76757 0.000808302
##
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
## state 1 state 2 state 3
## state 1 9.841498e-01 0.01585021 2.420180e-10
## state 2 1.471406e-44 0.98034212 1.965788e-02
## state 3 0.000000e+00 0.05615449 9.438455e-01
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3
## 9.999773e-01 9.870148e-06 1.283878e-05## Decoding state sequence... DONE
## Model AIC
## 1 sealHMMbd 3205.621
## 2 sealHMMb 3218.748
## 3 sealHMMdiurnDM 3235.313
## 4 sealHMMdiurn 3241.615## Computing pseudo-residuals...
While this combined model is better, given the complexity, I think I would stick to the bathymetry model for the moment and look into better ways to try to explain the data, especially since incorporating the time of day doesn’t fix the angle pseudo-residual problem.
This is where we stop the main tutorial. There are two other sections below. One that describes some of the mathematical details and one with activities that you can do on your own to further explore HMMs.
5 For the statistically/mathematically inclined
One of the most exciting aspect of the momentuHMM package is that you
can incorporate the potential effects of environmental covariates on the
probability of switching behavioural states. To understand how the
covariates can affect the behaviour, it’s useful to spend time
understanding the transition probabilities. This section is a bit more
technical, but I think is needed to fully understand what the estimated
covariate relationships mean.
If you have a model that has two behavioural states (state 1 and 2) and no covariates, you would have a transition probability matrix like this:
\[\begin{equation} \mathbf{\Gamma} = \begin{pmatrix} \gamma_{11} & \gamma_{12}\\ \gamma_{21} & \gamma_{22} \tag{5.1} \end{pmatrix} \end{equation}\],
where \(\gamma_{ij}\) represent the probability of switching from i to j. For example, \(\gamma_{11}\) represents the probability of staying in state 1 if you are already in state 1. Similarly, \(\gamma_{12}\) represents the probability of switching to state 2 if you are in state 1. By definition the row entries need to sum to 1. Which simply means that if you are in state 1 at time \(t-1\) you will end up in a state at time \(t\) (i.e., you don’t go in behavioural limbo).
Let’s look again at the transition probability matrix that was estimated in our seal example.
## Value of the maximum log-likelihood: -1810.643
##
##
## step parameters:
## ----------------
## state 1 state 2
## mean 0.05014675 4.902540
## sd 0.06428477 2.689235
##
## angle parameters:
## -----------------
## state 1 state 2
## mean 0.000000 0.000000
## concentration 2.789104 1.771595
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.346325 -4.840844
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2
## state 1 0.91264169 0.08735831
## state 2 0.00783846 0.99216154
##
## Initial distribution:
## ---------------------
## state 1 state 2
## 9.999995e-01 5.159521e-07We can see that the probability of switching in state 2 if you are in state 1 is 0.09 and so the probability of staying in state 1 is 0.91=0.91. The same logic applies to the probability of staying or switching from state 2. This means that for a model with \(N\) behavioural states, we only need to estimate \(N (N-1)\) transition probabilities (so 2 for the example here), because for each state at time \(t-1\) we will be able to calculate one of the transition probabilities based on the difference with the sum of the transition probabilities estimated (so here we can deduce \(\gamma_{11}\) based on \(\gamma_{12}\)).
To incorporate the covariates we are saying that the transition probabilities are changing through time as a function of the covariate values at that time step. Specifically,
\[\begin{equation} \label{eq:linkFx} \gamma^{(t)}_{ij} = \frac{\exp(\eta_{ijt})}{\sum^N_{k=1}\exp(\eta_{ikt})}, \tag{5.2} \end{equation}\]
where
\[\begin{equation} \label{eq:covRel} \eta_{ijt} =\begin{cases}\beta_{0ij} + \sum^p_{l=1} \beta_{lij} w_{lt} & \text{if} \ i \neq j,\\0, & \text{otherwise},\end{cases}\tag{5.3} \end{equation}\]
where \(w_{lt}\) is the \(l\)-th covariate at time \(t\) and \(p\) is the number of covariates considered. Equation is only a link function to insure that the transitions probabilities are between 0 and 1 and that the row sum to 1. The important relationship between the covariates and the transition probabilities are in equation , where the \(\beta_{0ij}\) is the base transition probability for the switch from state \(i\) to state \(j\) that is not affected by the covariates. The \(\beta_{lij}\) are the coefficient relating the transition probability to the \(l\)-th covariate. Note that \(\eta_{ii}\) (the probability of staying in state \(i\), the diagonal values in the transition probability matrix) are fixed to 0 because we only want to \(N (N-1)\) transition probabilities to be estimated and the row sums to one.
So in a model with no covariates, only the \(\beta_{0ij}\) for \(i \neq j\) (switching probabilities from state \(i\) to another state \(j\)) are estimated. If you look in the code box above you’ll see that we only have for 1 \(->\) 2 and 2 \(->\) 1, and that the transition probabilities for corresponding states are \(\gamma_{ij} = \frac{\exp(\eta_{ij})}{\sum^N_{k=1}\exp(\eta_{ik})}\).
For example, for the transition probability from state 1 to state 2.
## [1] 0.08735831# gamma_{12} calculated based on the equations above
exp(sealHMM2s$mle$beta[1])/(exp(0)+exp(sealHMM2s$mle$beta[1]))## [1] 0.08735831So the way we interpret the \(\beta\) values, is that if they are positive, it means that they increase the transition probability. Since the total value of \(\eta\) is what affects the transition probability, the magnitude and sign of the intercept (\(\beta_0\)) will affect when the transition probability goes from above or below 0.5.
This is a great way to incorporate covariates, but it means that it add \(N(N-1)\) parameters to your model, where \(N\) is the number of states. So if you have many behavioural states, it will results in a very complex model.
5.1 cosinor
The function cosinor internally creates the covariates
cos\((2 \pi \; \text{value}/\text{period})\) and
sin\((2 \pi \; \text{value}/\text{period})\). So for example
cos\((2 \pi \; \text{hour}/24)\).
6 Activities
6.1 Try it on your own data
Give it a try on your own data! If you don’t have data try the other activities.
6.2 Activity centers
This is fun but could be hard. Skip to the next activity if you find the material challenging.
Let’s look whether adding bathymetry improved the simulations. Here, it’s easy to have some problem, where the animal moves outside of the raster. So to help, we are setting the initial position as the initial position observed. We are also limiting the simulation to a time series of 200 steps. We are also setting the seed, because the animal will often wonder off map by chance.
set.seed(1)
sealSimB <- simData(model=sealHMMb, nbAnimals =1, spatialCovs = list(bathy=bathy),
initialPosition = as.matrix(sealreg[1,c("lon","lat")]),
obsPerAnimal = 100, retrySims = 10)##
Attempt 1 of 10...
Attempt 2 of 10...
Attempt 3 of 10...sealSimFitB <- fitHMM(sealSimB,
nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = par0b$Par, beta0=par0b$beta,
formula=bathyForm)## Warning in fitHMM.momentuHMMData(sealSimB, nbState = 3, dist = list(step = "gamma", : ginv of the hessian failed -- Error in svd(X): infinite or missing values in 'x'# Plot simulation
plot(water, col=bluePal(50))
points(sealSimB[,c("x","y")],
col=c("#E69F00", "#56B4E9", "#009E73")[viterbi(sealHMMb)],
pch=19, cex=0.5)
I’m not sure. The animal definitely wonders off in regions we wouldn’t expect it to go. But the state 1 is only in shallow waters.
Look at the momentuHMM vignette and see whether you can include
activity centers.
6.3 Modelling step length and turning angles
You can model the observations with a set of different distributions. To
change the distribution simply requires changing the argument dist in
the function fitHMM. In the example above, we used the gamma
distribution (using the name "gamma") for the step length. However,
you can use a variety of alternative distribution, including: Weibull
("weibull") and exponential ("exp"). Similarly, we used the von
Mises (named "vm") for the turning angle distribution, but one could
use the wrapped Cauchy ("wrpcauchy").
Explore whether changing the distribution affects the fit of the model.
Use AIC to compare the different models. Use the plot and viterbi
functions to see whether it has a meaningful biological effect.
For a description of the parameters needed for each of the
distributions, please see momentuHMM vignette.
#### Wrapped Cauchy for turning angle
kappa0wp <- c(0.9, 0.9) # Turning angle concentration parameter, kappa close
sealHMMwp <- fitHMM(sealPrep, nbState = 2,
dist=list(step="gamma", angle="wrpcauchy"),
Par0 = list(step=c(mu0,sigma0), angle=kappa0wp))## Model AIC
## 1 sealHMMwp 2878.146
## 2 sealHMM2s 3639.287#### Weibull for step length
mu0 <- c(0.6, 1) # Shape
sigma0 <- c(1, 5) # Rate
kappa0 <- c(0.9, 0.9) # Turning angle concentration parameter, kappa close
sealHMMwb <- fitHMM(sealPrep, nbState = 2, dist=list(step="weibull", angle="vm"),
Par0 = list(step=c(mu0,sigma0), angle=kappa0))## =======================================================================## Fitting HMM with 2 states and 2 data streams## -----------------------------------------------------------------------## step ~ weibull(shape=~1, scale=~1)## angle ~ vm(concentration=~1)##
## Transition probability matrix formula: ~1##
## Initial distribution formula: ~1## =======================================================================## DONE## Model AIC
## 1 sealHMMwb 3599.430
## 2 sealHMM2s 3639.287#plot(sealHMMwb) # Doesn't have a big effect
sum(abs(viterbi(sealHMM2s) - viterbi(sealHMMwb)))/nrow(sealPrep) # Only 27 (2%) different states## [1] 0.003577818#### Exponential for step length
mu0 <- c(2, 1) # Rate - has one less parameters
kappa0 <- c(0.9, 0.9) # Turning angle concentration parameter, kappa closesealHMMexp <- fitHMM(sealPrep, nbState = 2,
dist=list(step="exp", angle="vm"),
Par0 = list(step=c(mu0), angle=kappa0))## Model AIC
## 1 sealHMMwb 3599.430
## 2 sealHMM2s 3639.287
## 3 sealHMMexp 3923.112#plot(sealHMMexp) # It has a bit more of an effect
sum(abs(viterbi(sealHMM2s) - viterbi(sealHMMexp)))/nrow(sealPrep) # 128 different states (12%)## [1] 0.007155635Many of the distributions, including the gamma distribution, are
strictly positive. This means for example that you cannot have a step
length of exactly 0 if you simply use the gamma distribution. You can
however, use a zero-inflated distribution. Such a distribution assumes
that there is a probability \(z\) to observe a 0 and a probability \(1-z\)
to observe a positive value distributed according the the strictly
positive distribution of your choice (e.g., the gamma). This can be
achieved by adding a zero-mass starting value for each state in the step
parameters (see the help page of fitHMM under the argument Par0).
6.4 Effects of interpolation
The seal data we are using is not regular, by which I mean that the locations are not taken at regular time intervals. As a quick trick we regularised the data to 2 hours intervals using linear interpolation, but this is likely to affect the analysis. Here, explore the effects of using different time intervals. For example, look at the results when using 1, 4, 8 hours intervals. Look at both the fit of the model and the pseudo residuals. We commented out the plots for space, but do run them to visualise the effect of resolution on the data and results.
# Let's try different time intervals
# 1 hr
ti1 <- seq(seal$date[1], seal$date[length(seal$date)], by=60*60)
# 4 hr
ti4 <- seq(seal$date[1], seal$date[length(seal$date)], by=4*60*60)
ti8 <- seq(seal$date[1], seal$date[length(seal$date)], by=8*60*60)
# Interpolate the location at the times from the sequence
iLoc1 <- as.data.frame(cbind(lon=approx(seal$date, seal$lon,
xout = ti1)$y,
lat=approx(seal$date, seal$lat,
xout = ti1)$y))
iLoc4 <- as.data.frame(cbind(lon=approx(seal$date, seal$lon,
xout = ti4)$y,
lat=approx(seal$date, seal$lat,
xout = ti4)$y))
iLoc8 <- as.data.frame(cbind(lon=approx(seal$date, seal$lon,
xout = ti8)$y,
lat=approx(seal$date, seal$lat,
xout = ti8)$y))
# Create a new object that has the regular time and interpolated locations
sealreg1 <- cbind(date=ti1, iLoc1)
sealreg4 <- cbind(date=ti4, iLoc4)
sealreg8 <- cbind(date=ti8, iLoc8)
# # Quickly plot the regularized locations
# plot(seal$lon, seal$lat,
# pch=19, cex=0.5, xlab="Lon", ylab="Lat", las=1)
# points(sealreg1$lon, sealreg1$lat, pch=19, cex=0.5, col=rgb(1,0,0,0.5))
# points(sealreg$lon, sealreg$lat, pch=19, cex=0.5,
# col=rgb(0.8,0.8,0.8,0.5))
# points(sealreg4$lon, sealreg4$lat, pch=19, cex=0.5, col=rgb(0,0,1,0.5))
# points(sealreg8$lon, sealreg8$lat, pch=19, cex=0.5, col=rgb(0,1,1,0.5))
# Create the prep data
sealPrep1 <- prepData(sealreg1, type="LL", coordNames = c("lon", "lat"))
sealPrep4 <- prepData(sealreg4, type="LL", coordNames = c("lon", "lat"))
sealPrep8 <- prepData(sealreg8, type="LL", coordNames = c("lon", "lat"))
hist(sealPrep1$angle, breaks=40,
col=rgb(1,0,0,1), freq=FALSE,
ylim=c(0,2), las=1, border = NA)
hist(sealPrep$angle, breaks=40, add=TRUE,
col=rgb(0.8,0.8,0.8,0.8), freq=FALSE, border = NA)
hist(sealPrep4$angle, breaks=40, add=TRUE,
col=rgb(0,0,1,0.3), freq=FALSE, border = NA)
hist(sealPrep8$angle, breaks=40, add=TRUE,
col=rgb(0,1,1,0.3), freq=FALSE, border = NA)
## [1] 0## [1] 0## [1] 0## [1] 0# Ok let's fit the 3 state HMM again
# Fit a 3 state HMM
sealHMM3s1 <- fitHMM(sealPrep1,
nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu03s, sigma03s), angle=kappa03s))## =======================================================================## Fitting HMM with 3 states and 2 data streams## -----------------------------------------------------------------------## step ~ gamma(mean=~1, sd=~1)## angle ~ vm(concentration=~1)##
## Transition probability matrix formula: ~1##
## Initial distribution formula: ~1## =======================================================================## DONEsealHMM3s4 <- fitHMM(sealPrep4,
nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu03s, sigma03s), angle=kappa03s))## =======================================================================## Fitting HMM with 3 states and 2 data streams## -----------------------------------------------------------------------## step ~ gamma(mean=~1, sd=~1)## angle ~ vm(concentration=~1)##
## Transition probability matrix formula: ~1##
## Initial distribution formula: ~1## =======================================================================## DONEsealHMM3s8 <- fitHMM(sealPrep8,
nbState = 3,
dist=list(step="gamma", angle="vm"),
Par0 = list(step=c(mu03s, sigma03s), angle=kappa03s))## =======================================================================## Fitting HMM with 3 states and 2 data streams## -----------------------------------------------------------------------## step ~ gamma(mean=~1, sd=~1)## angle ~ vm(concentration=~1)##
## Transition probability matrix formula: ~1##
## Initial distribution formula: ~1## =======================================================================## DONE7 Thanks
Thanks to Théo Michelot for input and help!